88 


Journal of Glaciology, Vol. 57, No. 201, 2011 


Greenland ice sheet mass balance: distribution of increased mass 
loss with climate warming; 2003-07 versus 1992-2002 

H. Jay ZWALLY, Jun LI, Anita C. BRENNER, Matthew BECKLEY, Helen G. CORNEJO, 
John DiMARZIO, Mario B. GIOVINETTO, Thomas A. NEUMANN, John ROBBINS, 

Jack L. SABA, Donghui Yl, Weili WANG 

NASA Goddard Space Flight Center ; Code 614.1 , Green belt, Maryland 20771, USA 
E-mail: zwally@icesat2.gsfc. nasa.gov 

ABSTRACT. We derive mass changes of the Greenland ice sheet (GIS) for 2003-07 from ICESat laser 
altimetry and compare them with results for 1992-2002 from ERS radar and airborne laser altimetry. 

The GIS continued to grow inland and thin at the margins during 2003-07, but surface melting and 
accelerated flow significantly increased the marginal thinning compared with the 1990s. The net 
balance changed from a small loss of 7±3Gta -1 in the 1990s to 171±4Gta -1 for 2003-07, 
contributing 0.5 mm a -1 to recent global sea-level rise. We divide the derived mass changes into two 
components: (1) from changes in melting and ice dynamics and (2) from changes in precipitation and 
accumulation rate. We use our firn compaction model to calculate the elevation changes driven by 
changes in both temperature and accumulation rate and to calculate the appropriate density to convert 
the accumulation-driven changes to mass changes. Increased losses from melting and ice dynamics (17- 
206 Gta -1 ) are over seven times larger than increased gains from precipitation (10-35 Gta -1 ) during a 
warming period of ~2 K(10a) -1 over the GIS. Above 2000 m elevation, the rate of gain decreased from 
44 to 28 Gta -1 , while below 2000 m the rate of loss increased from 51 to 198 Gta -1 . Enhanced thinning 
below the equilibrium line on outlet glaciers indicates that increased melting has a significant impact on 
outlet glaciers, as well as accelerating ice flow. Increased thinning at higher elevations appears to be 
induced by dynamic coupling to thinning at the margins on decadal timescales. 


1. INTRODUCTION 

Mass changes in the Greenland ice sheet (GIS) are of 
considerable interest because of their sensitivity to climate 
change, the mass-loss contribution to sea-level rise, and the 
likelihood of an increasing rate of mass loss with climate 
warming. In recent years, the GIS has experienced increased 
surface melting from increasing temperatures (Hall and 
others, 2008; Box and others, 2009), increased snow accu- 
mulation from increasing precipitation (Box and others, 
2006; Hanna and others, 2008), accelerated ice flow from 
melt/acceleration effects (Zwally and others, 2002b; Joughin 
and others, 2008b; Van de Wal and others, 2008; 
Bartholomew and others, 2010) and increased discharge 
from accelerating outlet glaciers due to warming ocean 
waters (Holland and others, 2008; Straneo and others, 
201 0). In response, the GIS has been thinning at the margins 
and growing in its interior (Krabill and others, 2000; 
Johannessen and others, 2005; Zwally and others, 2005; 
Pritchard and others, 2009). The net rate of ice loss appears 
to have accelerated in recent years compared with the 
1990s, but published values of the rate of mass loss and 
changes with time have differed widely (Alley and others, 
2007; Shepherd and Wingham, 2007; AMAP, 2009). 

Estimates of input and output (AMAP, 2009) for the GIS 
for a state of near balance give a mass gain of 550 Gta -1 
from precipitation (less evaporation), which is balanced by 
mass loss through meltwater and blowing snow (370 Gta -1 ) 
and glacier discharge (180 Gta -1 ). However, these values 
have large uncertainties and significant interannual vari- 
ations. Variations in precipitation, temperature and other 
factors cause variations in the mass exchange at the surface 


and in the dynamics of ice flow toward the margins on land 
and through outlet glaciers into the ocean. 

A review (Shepherd and Wingham, 2007) said 'it is 
reasonable ... to conclude the GIS is losing about 
100 Gta -1 ', which was consistent with the high-resolution 
analysis of Gravity Recovery and Climate Experiment 
(GRACE) data (Luthcke and others, 2006) that gave a net 
mass loss of 101 Gta -1 for the period 2003-05. Furthermore, 
a 122 ±47 Gta -1 increase in loss rate by 2005 compared 
with 1996 based on a mass flux analysis (Rignot and 
Kanagaratnam, 2006) is approximately consistent with the 
increase shown by the GRACE results for 2003-05 
compared by Luthcke and others (2006) with the near 
balance for 1992-2002 from satellite-radar and airborne- 
laser altimetry (Zwally and others, 2005). 

This paper focuses on the GIS mass balance for 2003-07 
and the changes since 1992-2002, including the spatial 
distributions of mass balance and changes by drainage 
system and elevation. Our 2003-07 results are derived from 
the Ice, Cloud and land Elevation Satellite (ICESat)'s laser 
altimeter measurements (Zwally and others, 2002a) of 
surface elevations using some of the methodology previously 
applied (Zwally and others, 2005) to data from European 
Remote-sensing Satellite (ERS) radar altimetry and airborne 
laser altimetry (Airborne Topographic Mapper (ATM)). Sig- 
nificant improvements for determining mass changes, dM/dt, 
from elevation changes, dH/dt, for both the 1992-2002 and 
2003-07 periods are (1) calculation of the appropriate 
densities to apply to changes in the thickness of the firn-ice 
column and (2) partitioning of dH/dt and dM/dt into 
components driven by accumulation variations and those 
driven by melting and ice dynamics. For 1 992-2002, we use 


Zwally and others: Greenland ice sheet mass balance 


89 



Fig. 1. Schematic of elevations h(Xj, a, tj) for successive data 
collected on repeat tracks at times t-, along an ICESat reference track 
with cross-track slope a and distance x,- from the reference track, h 
is shown increasing with time. 


the previous dH/dt (Zwally and others, 2005) calculated from 
elevation changes, dh/dt, calculated at crossovers where ERS 
satellite tracks cross on consecutive ascending and descend- 
ing passes, with ATM data (Abdalati and others, 2001) added 
to fill in gridcells on outlet glaciers. For 2003-07, we 
calculate dH/dt from dh/dt calculated from ICESat data 
collected on successive satellite passes along closely spaced 
ground tracks ('repeat tracks'), which significantly increases 
the along-track spatial resolution and the density of derived 
dh/dt compared with crossover analysis. 

2. ELEVATION CHANGE, dh/dt, FROM REPEAT- 
TRACK ANALYSIS OF ICESAT DATA 

Our repeat-track analysis to derive dh/dt from ICESat data is 
enabled by ICESat's capability to point the laser beam to 
reference tracks over the ice sheets to ±100m (1 cr). 
However, because the tracks do not repeat exactly, the 
surface elevation, h(x, y, tj ), measured along a repeat pass (y 
direction) at time tj depends on the cross-track surface slope, 
a, and the cross-track distance, x, as well as actual elevation 
changes with time. This mixing of apparent height changes 
caused by cross-track displacements and real height 
changes, h(t,), is illustrated in Figure 1 and by 

h(x, y n a, tj) = xtan a + t, (^j + h 0 (y r , to), (1 ) 
where h 0 is the elevation at the position, y r , of the reference 


track at t 0 . Use of Equation (1) assumes the cross-track 
surface slope is constant over distances of ~200 m and that 
dh/dt is constant with time, i.e. h(t) changes are linear. 
Although important seasonal and interannual variations in 
h{x, y a, t) are ignored in this solution, departures from 
linearity are nevertheless retained in the height residuals 
relative to the linear solution and can be further analyzed. 
Measured h{x, y a , t) are first interpolated to equally spaced 
(172 m) reference points along each track. We then solve 
Equation (1) by least-squares methods for three parameters, 
a , dh/dt, h 0 , at each reference point, which requires data on 
at least N=4 along-track repeats. For solutions with N=10 
repeats, for example, we obtain an accuracy, cr^h/dt, for dh/dt 
of better than 25 cm a -1 for a <2° at one-point along-track 
solutions (Fig. 2a). The resulting standard errors, cr dH / df and 
cr a , on dh/dt and a decrease with increasing Nby 1/(N-3) 1/2 
as expected. We use data from 12 laser campaigns (also 
designated as Laser 2a, 2b, 2c, 3a, 3b, 3c, 3d, 3e, 3f, 3g, 3h 
and 3i) from the fall of 2003 to the fall of 2007, i.e. F03, 
W04, S04, F04, W05, S05, F05, W06, S06, F06, W07 and 
F07, where F indicates boreal fall, W indicates boreal 
winter, S indicates boreal spring and 03 indicates 2003, etc. 
The effective length of each campaign was 33 to ~36days. 
(F03 was 45 days, but the additional ground tracks were not 
covered in subsequent campaigns.) However, data are not 
continuous along each track due to losses in thick-cloud 
cover, so N varies with location from a maximum of 12. The 
data versions are all release 428, except for S04, which is 
release 128. Data can be obtained from NASA's archive at 
the US National Snow and Ice Data Center (http://nsidc.org/ 
data/icesat/i ndex.html). 

The dh/dt results are improved by a seven-point solution 
that uses data at three points on either side of each 
reference point, i.e. strips that are ~1 km along track, and a 
quadratic surface in the along-track y direction that is 
derived from multiple repeats. The maximum number of 
h(x, y, a, t) values used in the seven-point solutions is 84. 
We obtain seven-point solutions at 97% of the locations. 
For the other 3%, which are mostly near the edge of the ice 
sheet where a satisfactory seven-point solution cannot be 
obtained, we use one-point solutions. The formal cr d/7/df on 
dh/dt from the seven-point solutions are multiplied by \fl 
to account for the non-independence of consecutive along- 
track solutions. 
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Fig. 2. (a) Error on ICESat dh/dt (ma ] ) as a function of the number of repeats, N, used in the solution at each along-track location separated 
by 1 72 m. 0.25 m a -1 accuracy is achieved for N> 10 repeats for slopes <2°. Dashed curves are 1/(/V-3) 172 fits to the errors, (b) Error on 
slope a as a function of N. Slope error is about ±20% for N= 10. 
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Fig. 3. (a) ICESat surface elevation {h) profiles (ten-point smoothing) 
relative to the first profile from fall 2003 across Jakobshavn Isbras 
along T419 (Fig. 14). Maximum surface lowering is 27.6m in 
4 years. F03 is fall 2003, W04 is winter 2004, S04 is spring 2004, 
etc. (b) Profiles of h from 1030 to 1465 m, derived dh/dt, and crdh/dt 
from along-track solutions at 172 m spacing. Units of h are from 
1000 m (0 on y-axis) to 1500 m (5 on y-axis). Maximum dh/dt is 
7 m a -1 near the center of the glacier. Matching peaks in h and dh/dt 
indicate oscillations are growing with time, (c) h(t) at crossings of 
ICESat tracks with a flowline near the center of the glacier (Fig. 1 4). 
Over 4 years, h(t) is close to linear. 

We further improve the relative accuracy to the cm level 
for each of the 12 campaigns by referencing to the ocean 
surface, i.e. the open water within sea ice, on the Arctic 
Ocean. This reduces inter-campaign range biases. We 
empirically determine an elevation correction, d, from the 
time variation of the average mean sea surface over the 
Arctic Ocean using the methods described by Zwally and 
others (2008). Although the magnitude of d is spatially 
variable, its time variation appears to be the same over the 
Arctic Ocean, so we take it to be applicable to Greenland as 
well. The respective values of d for the 12 periods are 


-0.049, 0.045, 0.132, -0.061, -0.104, 0.021, 0.011, 0.014, 
0.012, 0.074, 0.004 and 0.043 m. We reduce the time 
variation of these d values by 0.003 m a -1 to account for the 
current rate of sea-level rise, and then subtract the reduced d 
values from the measured elevations. The linear trend in the 
reduced d is 0.006 ma _1 , which averaged over all of 
Greenland increases the overall mass loss by ^9Gta -1 
compared with data without the d correction applied. 

2.2. Characteristics of dh/dt on Jakobshavn Isbrae 

We select Jakobshavn Isbrae, west-central Greenland, to 
show the characteristics of the dh/dt derived from ICESat- 
measured surface elevations, /?, in Figure 3a, including the 
small cr d /yd t relative to dh/dt, the details of the oscillations in 
the h and dh/dt profiles, the consistency of the oscillations 
on consecutive profiles, and the small variation from 
linearity in the h(t) at different locations. Jakobshavn is of 
particular interest because of its rapid speed-up and 
increased thinning rate inland of the grounding line as the 
floating ice tongue thinned and broke away during 1996- 
2004 (Joughin and others, 2008a; Thomas and others, 2009). 
Profiles of h for 12 ICESat repeats along track 419 (T419) 
relative to the first pass in fall 2003 are shown in Figure 3a 
along with dh/dt, its formal error, crdh/dt, and elevation in 
Figure 3b (track locations are shown in Fig. 14 further 
below). 

The largest dh/dt derived in this region is 1 1 .3 m a -1 at the 
end (69.1 7° N, 49.32° W) of the discontinuous T85, near the 
center of the glacier downstream from T419. The linear 
character of h(t) in Figure 3c at six crossings of ICESat tracks 
with a flowline down the center of the glacier (Fig. 14), 
where dh/dt values range from -8.7 to -0.1 ma _1 below and 
above the EL, shows that the glacier's thinning rate is nearly 
constant during 2003-07. Therefore, any residual effects 
from removal of the ice tongue do not appear to be causing 
either a significant increase or decrease in the thinning rate 
in this portion of the glacier. 

3. MAPPING OF dH/dt 

We average the derived ICESat dh/dt values over ~50km 
gridcells to obtain the average elevation change, dH/dt, in 
each cell (Fig. 4a), along with the formal error, (Jd H /dt, of the 
average dH/dt (Fig. 5a). The dH/dt values in Figure 4a are 
similar to the map of along-track dh/dt derived from ICESat 
data as shown in figure 2 of Pritchard and others (2009), 
which described the patterns of dynamic thinning but did 
not calculate dM/dt, as we do in section 4. The ICESat (TdH/dt 
are approximately one-third of the corresponding values of 
those for ERS radar altimetry (Fig. 5b), which are computed 
from the time-series analysis of crossover differences (Zwally 
and others, 2005). The smaller cr dH/df for ICESat reflects the 
higher accuracy of the laser versus radar altimetry and 
possibly the higher accuracy of the repeat-track analysis 
(even though data collection was not continuous). 

The dM/dt, area and other parameter calculations in the 
following sections are made using the 50km gridcell 
averages. However, at the edges of the ice sheet, the 
50 km gridcells are split into portions lying on and off the ice 
sheet, using an ice-sheet boundary with 1-2 km resolution. 
For the dH/dt from ICESat data, only dh/dt data within the 
portion of the cells lying on the ice sheet are used in the 
calculations. For the dH/dt from ERS/ATM, the procedures 
for selection of crossovers, the use of ATM data, and optimal 
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Fig. 4. (a) ICESat dH/dt= (dh/df) averaged in ~50km gridcells. (b) Temperature-driven part of the firn compaction dC r /df times 10. (c) d // 
dt=(dH/dt)-(dC 1 /dt)-{dB/dt). Although dC-j/dt is generally smaller than dH/dt, it shifts the dH/df=0 line by ~100km in the northwest as 
shown by the arrows. In addition to the delineation of the drainage systems (black curves), the 2000 m contour (dotted curve), the EL (dashed 
curve) and the ice-sheet boundary (red curve) are shown here and in Figures 5-9. 


interpolation to fill in missing cells were given by Zwally 
and others (2005); here we also weight all the boundary 
cells by the fraction of the cell area lying within the 
boundary. An important correction to dH/dt from radar 
altimetry is the backscatter correction previously applied 
(Zwally and others, 2005), which empirically corrects for 
variations in the effective backscattering depth in the firn 
and is similar to that applied by Davis and others (2005) and 
Wingham and others (2006). 


4. DERIVATION OF ICE-MASS CHANGES, dM/df, 
FROM ELEVATION CHANGES, dH/dt 

Derivation of dM/dt from dH/dt requires consideration 
(Zwally and Li, 2002) of the surface mass balance, firn 
compaction, dynamics, and underlying bedrock motion: 


dH(f) 

df 


m 

P sf 


- v fc (0 - 


A )(0 

Pi 


-v + **, 

df ' 


(2) 


where A(t) is the snow accumulation rate, p s f is the density of 
near-surface firn (typically 0.33), Vf c (f) is the velocity of firn 
compaction at the surface, A b (f) is the ablation rate, p\ = 0.90 
is the density of glacier ice, V lce is the vertical velocity of the 
ice at the firn/ice transition, and dB/dt is the vertical bedrock 
motion. In steady state with a long-term average accumu- 
lation rate of (A), Vj ce = (A)/p\ in the accumulation zone and 
Vf c = 1 .7 \/i ce using p sf = 0.33 and pj = 0.90. The units of A(t) 
and A b (f) are cmw.e. a -1 , the units of velocities and 
elevation changes are cm a -1 , and all densities are relative 
to Pwater = 1 . H, A and dB/dt are positive upward and Vf c , A b 
and Vj ce are positive downward as commonly defined for 
these parameters. 

For non-steady state, we rewrite the ice-sheet terms in 
Equation (2) as perturbations from steady state, assuming 


that dB/dt is constant during the time of the measurements: 

dH dH a d Cat dH b dH d dB 

df df df df df df ' K J 

where dH a /df=(A(f)-(A))/p s f is the direct change caused by 
A(f), dC>\ 7 /df=-(\/f c (f)-(Vf c )) is the change in firn-compac- 
tion rate driven by both A(f) and temperature T(f), dH b /dt=- 
(A b (t)-(A b ))/p\ is driven by changes in the ablation rate, dH d / 
df=-(Vj Ce (f)-(Vj ce )) is driven by dynamic changes in the ice 
flow relative to the longterm, (A), and () indicates long-term 
averages of the various parameters. All terms in Equation (3) 
are positive upward and are defined as averages over some 
area, which in our following analyses are 50 km gridcells. 

Positive values of dH d /dt indicate that the ice flow is less 
than that required to balance the long-term average (A), 
which is described as dynamic thickening in the same 
manner as negative dH d /dt is described as dynamic thinning. 
EHowever, dH d /dt includes not only long-term changes in the 
ice dynamics, but also short-term changes. For example, if 
dH d /dt is determined during two different periods at an 
interval of 10 years, then both values of dH d /dt would 
include multi-decadal dynamic imbalances compared with 
the multi-decadal average (A), and their difference would 
indicate the sub-decadal dynamic changes (see section 6.3). 

4.1. Separation of accumulation and temperature 
effects on firn compaction 

Thefirn-compaction (dC A7 /dt) and dB/dt terms in Equation (3) 
change dH/dt but not dM/dt. Both dH a /dt and dC A1 /dt are 
affected by A(f), and dC A7 /dt depends strongly on 7(f) as 
previously used (Zwally and others, 2005) to calculate dM/ 
df. Importantly, dC A1 /dt depends on the time history of both 
A(f) and 7(f) as their effects propagate into the firn. We 
assume the compaction effects of A(f) and 7(f) are separable 
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Fig. 5. (a) Maps of cr dH / d? for ICESat 2003-07 computed from the cell-average dH/df of the dh/dt solutions at the individual reference points, 
(b) (JdH/dt for ERS 1992-2002 is computed from time-series analysis of crossover differences (Zwally and others, 2005). 


so that 


dC AT _ d Ca dC T , . 

df ~ dt df U 

which should be valid at least for small A(t) and 7(f) 
perturbations. The total elevation change from both A(t) and 
7(f) perturbations becomes 


dH a cAT _ d/~/ a d Ca dC 7 = dH a CA dC 7 
dt ~ dt dt dt dt df 


where dH a CA /dt is the total change caused by A{t), 
including the direct change and the change induced in the 
compaction rate. Equation (3) then gives 

dl dH dC T dB dH a CA dH b dH d 

S s dT~-d7 = — + 1F + 1F' (6) 


where d//df is an approximation of the change in thickness 
of the ice column, defined by treating dB/dt and dCj/dt as 
corrections to dH/dt. 

In Equation (6), the three terms on the right-hand side 
involve changes in mass, and 


dM 

dt 


Pa 


dH a C A 

dt 


+ p\ 


dH bd 

dt 


Ar f 


Pa dH a CA , (dl dH a CA \ 
df +p \dt d T ) 


Ar f, 


(7) 


where p a is the appropriate density for the A(f)-driven 
changes, p\ is the density of ice appropriate for the ablation 
and dynamic terms (combined as dH bd /dt= (dH b /dt) {dH d / 
df)), Ar is the area of a firn/ice column or map gridcell, and 
f is the fraction of the cell inside the ice-sheet boundary. We 
use pi = 0.90 for the typical density of glacier ice in the 
ablation zone and 0.91 in the accumulation zone, taking 
into account small variations from the 0.91 7 density of solid 
ice due to air content. dM/dt is also used in units of mass 
change per unit area, which is defined as the surface flux 
equivalent (SFE). The SFE is equal to the negative of the 
mass-flux change at the surface that would be required 
for the ice column to be in balance. Maps of dM/dt 


are presented in SFE units of Gta 1 (1 00 km) 2 , where 
Gta -1 (100 km) -2 = 100 kg a -1 m“ 2 . 

4.2. Separation of accumulation and melting/ice- 
dynamic driven mass changes 

An important aspect of our formulation is the separation of 
total dM/dt into two components: dM a /dt= p a dH a c A /dt, 
which is driven by A(t) perturbations, and dM bd /dt= p\(dH bd / 
df), which is driven by ice dynamics and by ice ablation 
below the equilibrium line (EL). The EL varies from ~1 000 m 
in the south, to 1 600 m around 72° N, to 200 m in the north 
(Zwally and Giovinetto, 2001 ). Above the EL, dhl b /dt= 0 and 
we use d H bd /d t = d H d /d f. 

The primary value in Equations (6) and (7) is dH/dt 
measured over some time, Af. For dB/dt, we use the radial 
components of three models of isostatic rebound (Ivins and 
others, 2001; Huybrechts, 2002; Peltier, 2004), labeled (dB/ 
df)| H/ p respectively, interpolated to our gridpoints. As in 
Zwally and others (2005), we weight these models to 
account for distributions that showed similar patterns, 
specifically dB/dt= [(dB/dt) P /4] + [(dB/dt) H / 4] + [(dB/dt)\/ 2]. 
We obtain the values of p a , dH a cA/dt and dC T /dt from our 
firn-compaction model (Zwally and Li, 2002), which was 
modified to include the effects of surface melting and 
percolation (Li and others, 2007) and is further modified 
here to include a variable accumulation rate, A(t), as 
described in more detail by Li and Zwally (in press). The 
model is driven by a T(t) temperature record from 1982 to 
the present as derived from Advanced Very High Resolution 
Radiometer (AVHRR) satellite measurements (Comiso, 2003) 
with values extended and updated using new calibrations. 

We parameterize the accumulation perturbations, 
A(t)-(A), as a function of 7(f) using a sensitivity of 
(5 =L 1 .5)%8A K -1 . The Intergovernmental Panel on Climate 
Change (IPCC) report (Solomon and others, 2007) summar- 
ized the sensitivity of accumulation rate to temperature with 
a range of 4. 7-8.5% 5A K _1 . Data from Greenland ice cores 
(Clausen and others, 1988) give a range of 4-6% 5A KT 1 using 
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Fig. 6. (a) Distribution of d//df for 1992-2002 derived from ERS/ATM (Zwally and others, 2005). (b) Distribution of d//df for 2003-07 from 
ICESat. (c) Change in dl/dt between 1992-2002 and 2003-07. Increased thinning is evident around most of the margins. Increased 
thickening is shown in the northern drainage systems (DS1 and DS2) and at higher elevations in the south (DS5 and DS6). The central region 
has a small thinning above ~2000m in the latter period. 


a sensitivity of 5 1s O to temperature of 0.69%o K -1 (Zwally and 
Giovinetto, 1997). For 1988-2005, our T{t) data and A(T(t)) 
parameterization give an A(f) increase of 0.6% a -1 . This rate 
is close to the 0.7% a -1 we calculate for 1988-2004 from 
results of a regional-climate and surface mass-balance model 
using linear fits to the trends in the percolation and dry-snow 
zones given in figure 1 1 a of Box and others (2006) with zone 
areas weighted by 30% and 70%. 

4.3. Density, p a , for accumulation-driven mass 
changes 

Our compaction model first calculates dH a cAr/dt using 
both T(t) and A(t) f and then calculates dCj/dt using T(t) with 
constant (A). Equation (5) then gives the d/-/ a o\/df, and 
dHbd/df is the residual part of dH/dt according to Equa- 
tion (6). Equation (6) gives dl/dt for calculation of dM/dt in 
Equation (7), with all values calculated as averages over At. 
The average p a is calculated over Affrom p a = AM A /AH a cA, 
where AM a and A H a cA are the integrals over At of 
(A(f)-(A)) and dH a o\/df. Our calculated p a is 0.525 for 
1992-2002 and 0.527 for 2003-07, averaged over the ice 
sheet above the EL. The distribution of p a in Figure 7c shows 
that p a is generally smaller (e.g. <0.5) in the northern part of 
the ice sheet and toward higher elevations where lower A 
reduces the rate of firn compaction. The largest values of p a 
(e.g. >0.65) are in the areas with high A at lower elevations 
mostly in the southern part. 

Previous estimates of dA4/df from dH/dt made by various 
authors assumed there is an effective density, p e ff, so that 
dA4/df=p e ffd//df, or when firn compaction terms are 
neglected dM/dt= p e ^(dH/dt-dB/dt). For example, we pre- 
viously used a density of 0.90, which is appropriate for solid 


ice if all the dH/dt is from long-term changes, and we 
illustrated a lower dM/dt estimate using p = 0.4, which is 
appropriate if all the dH/dt is from the addition or depletion 
of the upper part of the firn (Zwally and others, 2005). Other 
authors used values ranging from 0.35 to 0.91 7 (e.g. 0.35 in 
Davis and others (2005); 0.300 for the Interior' and 0.900 to 
'seaward' in Thomas and others (2006); and 0.350 and 
0.91 7 in Wingham and others (2006)). EHowever, Equation (7) 
shows that p e ff is strictly valid only if A(f) is constant, so that 
dH a cA/dt = 0 and p e ff = 0.9. Calculations of dM/dt using 
Equation (7) show that in some locations p e ff values 
calculated using (dM/dt)/(dl/dt) are unrealistically outside 
the 0.3-0.91 range for firn/ice densities. Therefore, using a 
fixed value of p e ff is only valid if A(f) variations are 
insignificant (i.e when dH a cA/dt = 0 and p e ff = 0.9). Consist- 
ent with this conclusion, EHelsen and others (2008) used a 
compaction model based on the model of Zwally and Li 
(2002) and showed that a significant portion of the dH/dt 
observed in East Antarctica can be caused by temporal 
variations in accumulation. 


5. RESULTS 

The spatial distributions of our calculated dCj/dt and dl/dt 
for 2003-07 are given in Figure 4b and c. Although dCj/dt is 
generally smaller than d H/dt, the larger values of dl/dt 
evident at higher elevations, in particular, are a result of 
warmer temperatures increasing the rate of firn compaction. 
We also recalculate dCj/dt and d//df for the 1992-2002 
period using the dH/dt from Zwally and others (2005), and 
compare the d//df for the two periods in Figure 6a and b, 
with their difference in Figure 6c. To compare the rates of 
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Fig. 7. Maps for the 2003-07 period, (a) Accumulation-driven elevation change, dH a cA/dt. (b) Ablation- and dynamic-driven elevation 
change, dH bd /dt. (c) Relative density, p A , of the firn for the dH a cA/dt component. 
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Fig. 8. Map of dM/dt (Gta -1 (100 km) -2 ) for 2003-07 derived from 
ICESat data. The numbers around the map are the net gains (red) or 
losses (blue) (Gta -1 ) for the DS1-DS7. The numbers in parentheses 
are the changes from the 1992-2002 period. All drainage systems 
except DS5 are divided into subsystems (e.g. DS1.1, etc.). The 
largest increases in the rates of mass loss are in the west central 
region and the southeast. The northern DS1 and DS2 are gaining 
mass. Most of the areas above and below 2000 m (dotted curve) are 
gaining and losing mass, respectively. Mass losses below the EL 
(dashed curve) are from a combination of increased melting and 
acceleration of outlet glaciers. 


ice-sheet thickness change, the comparison of d//df is more 
meaningful than comparison of dH/dt maps, because of the 
magnitude of dCj/dt that changed with climate warming 
between the respective periods. 

During 1992-2002, much of the ice sheet was shown to 
be thinning below 2000 m elevation and thickening above 
2000 m. During 2003-07, the largest changes in d//df are the 
increases in thinning along most of the ice-sheet margin 
below the EL and below 2000 m. Above 2000 m, increased 
thickening is observed in the southwest drainage systems 
(DS5 and DS6) and in the northeast (DS1 .2, DS1 .3 and DS2). 
EHowever, this increased thickening is mostly counterba- 
lanced by a large area of small thinning in the central region. 

The distributions of d/-/ a o\/df, dH bd /dt and p a for 2003- 
07 are shown in Figure 7a-c. The resulting distribution of 
dM/dt is in Figure 8, and the map of the calculated error, 
crdM/dt, is in Figure 9a, along with the cr dM / dt for 1992-2002 
in Figure 9b. cr dAd/df is calculated using propagation of errors, 
<jy - £(<9y/<9x/cr/) 2 , which applied to Equation (7) gives 

(^d/u/df) - (Arfpi) 2 |pi 2 [o- dH / df 2 + crdcr/df 2 + ^de/df 2 ] 

dH a CA 

( 8 ) 

For ICESat, cr dH /df in Figure 5a is computed from the cell- 
average dhl/dt of the dh/dt solutions at the individual 
reference points. For ERS, cr dH / dt in Figure 5b is computed 
from the time-series analysis of crossover differences (Zwally 
and others, 2005). cr dc T / dt is calculated from the fit of Cj{t) 
for the appropriate period (October 2003-October 2007 for 
ICESat and April 1992-October 2001 for ERS), as the un- 
certainty in the time derivative in Cjit) = a + (dC T /dt)t+ 
csin(caf) + d cosfcaf), where cd = 27r/365.25 and Ms in days. 
o'dB/dt is computed as half the difference between the similar 
and different values used to compute dB/dt. (JdH*CA/dt is the 
difference between d/-/ a o\/df for the cases where 


+ [(Pi — Pa)o'dHaCA/dt] + 
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Fig. 9. (a) Maps of <JdM/dt (a) for 2003-07 from ICESat data and (b) for 1992-2002 from ERS and ATM data. 


accumulation changes at a rate of 5% K -1 and 8% K 1 . cr pa is 
the standard deviation of the distribution of p a . 

The values of dH/dt, dCj/dt, dl/dt, dH a C A/dt and dH bd /dt 
as functions of elevation averaged over the ice sheet are 
shown in Figure 1 0. The importance of dCj/dt as a correction 
to dH/dt is illustrated by the close agreement between dl/dt 
values (5.2 vs 4.9 cm a -1 ) for 1992-2002 and 2003-07 
respectively at elevations above 2200 m, in contrast to the 
large difference between the dH/dt values (4.9 vs 1 .0 cm a -1 ). 
Without dCy/df, which is larger for the warmer 2003-07 
period, the dH/dt alone would have incorrectly suggested 
that the ice-sheet thickening of ^5 cm a" 1 at higher eleva- 
tions (Johannessen and others, 2005; Zwally and others, 


2005) had significantly diminished. The average dCj/dt over 
the accumulation zone are -0.2 and -4.1 cm a -1 for 1992- 
2002 and 2003-07 respectively, and the associated volume 
changes of -3 km 3 a -1 and -59 km 3 a -1 would lower the 
derived dM/dt by 3 Gta -1 and 54 Gta -1 if the dCj/dt 
correction to dH/dt were not applied. 

The average dB/dt is 0.06 cm a -1 , which adjusts dM/dt by 
-1 Gta -1 . The close agreement of dl/dt at high elevations 
averaged over the ice sheet between the two periods 
suggests that the laser and radar results are actually quite 
comparable, despite the time difference between the 
measurements and technical differences between the laser 
and radar measurements. There is also close agreement 


Table 1. dM/dt (Gta ^ by drainage system (DS): four DS in north and southwest (N SW; 1, 2, 5, 6), four DS in west central and southeast (WC 
SE; 3, 4, 7, 8) and totals. (Totals are rounded to integers after summation of values by DS before rounding) 


DS 

8 

7 

6 

5 

DS 

4 

3 

2 

1 

N SW 

WC SE 

Total 

2003-07 

Total 

-33 

-14 

-4 

-10 

-75 

-51 

13 

3 

3 

-174 

-171 ±4 

>2000 

5 

8 

11 

1 

-10 

-6 

13 

5 

31 

-3 

28 ± 1 

<2000 

-38 

-23 

-15 

-10 

-65 

-45 

0 

-2 

-28 

-171 

-1 98 ±4 

>EL 

-9 

0 

8 

-3 

-39 

-25 

17 

5 

27 

-72 

-46 ±2 

<EL 

-24 

-15 

-12 

-6 

-37 

-26 

-4 

-2 

-24 

-101 

-125 ±3 

1992-2002 

Total 

3 

6 

1 

-3 

-9 

-9 

6 

-2 

1 

-8 

-7 ±3 

>2000 

11 

12 

6 

0 

0 

3 

9 

4 

18 

26 

44 ±1 

<2000 

-7 

-6 

-5 

-3 

-9 

-12 

-3 

-6 

-16 

-35 

-51 ±3 

>EL 

9 

11 

6 

-1 

-3 

-2 

7 

0 

12 

15 

26 ± 3 

<EL 

-6 

-6 

-5 

-2 

-6 

-6 

-2 

-2 

-10 

-23 

-33 ±1 

Change from 1992-2002 to 2003-07 
Total -36 -20 

-5 

-6 

-66 

-42 

7 

5 

2 

-165 

-164 ±5 

>2000 

-6 

-4 

5 

1 

-10 

-9 

5 

1 

13 

-29 

-16±1 

<2000 

-30 

-16 

-10 

-7 

-56 

-33 

2 

4 

-11 

-136 

-147 ± 5 

>EL 

-18 

-11 

3 

-2 

-35 

-22 

10 

5 

15 

-87 

-72 ±3 

<EL 

-18 

-9 

-7 

-4 

-31 

-20 

-2 

0 

-14 

-78 

-92 ±3 
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Fig. 10. All parameters averaged in 500 m elevation bands over the 
ice sheet versus elevations for 1992-2002 and 2003-07. Averages 
are over all elevations weighted by area, (a) Measured dH/dt 
decreased with time at all elevations, (b) Firn-compaction correc- 
tion dCy/df driven by temperature variations only, (c) dl/dt={dH/ 
df) - (dCy/df) - (dB/dt) shows little change with time above 2200 m. 
(d) dH a cA/dt, the portion of d//df driven by accumulation 
variations, increased slightly at all elevations, (e) dH bd /dt, the 
portion of d//df driven by changes in ablation (below EL only) and 
ice dynamics, decreased at all elevations. The increased thinning 
below the EL is from increased melting and acceleration of outlet 
glaciers. The increase in dynamic thinning extends inland to the 
highest elevations. 


(Fig. 11) between the two periods in both the derived 
dA/fbd/df and the dA 4/dt at all elevations in four (DS1, DS2, 
DS5 and DS6) of the eight drainage systems, with average 
values being slightly higher in DS1 and DS2 and slightly 
lower in DS5 and DS6 in the latter period. 

The distribution of dM/dt (Figs 8 and 1 1 ; Table 1 ) shows 
that during 2003-07 most of the GIS above 2000 m was 
gaining mass (+28 ±1 Gta -1 ) and most of it below 2000 m 
was losing mass (-1 98 ±4 Gta -1 ). Both these rates are 
lower than the 1992-2002 values of +44 ±1 Gta -1 and 



Fig 11. dA 4/dt is total rate of mass change, dMJdt is the component 
driven by temporal variations in snow accumulation, and dAd^d/dfis 
the component driven by ablation and ice dynamics, all averaged 
by 500 m elevation bands over the ice sheet for the 1 992-2002 and 
2003-07 periods. Circled symbols are totals for all elevations 
weighted by area. 


-51 ±3 Gta -1 . The total loss of 171 ±4Gta -1 during 2003- 
07 is in close agreement with two GRACE-based loss rates 
for similar time periods of 1 58 =t 8 Gta -1 for October 2003- 
November 2007 (Luthcke and others, 2009) and 1 79 =L 
25 Gta -1 for February 2003-January 2008 (Wouters and 
others, 2008). Our derived mass loss of 171 ±4 Gta -1 
represents an increased loss rate of 1 64 ± 5 Gta -1 compared 
with the near-balance conditions during 1992-2002, for 
which we obtain a revised small loss of 7±3Gta -1 . The 
previous small gain (Zwally and others, 2005) of 1 1 Gta -1 is 
revised here due to several factors: -13.8 Gta -1 change is 
from improvements in the firn-compaction model, 
-9.2 Gta -1 from use of recalibrated T(t) data, +11. 9 Gta -1 
from changes in the grid-map projection and use of partial 
gridcells at the ice edge, and 7.8 Gta -1 from partitioning of 
the height changes into accumulation-driven changes with 
an average density of 0.53 and ablation-/dynamic-driven 
changes with an average density of 0.9. 

Our partitioning of the dM/dt gives an accumulation- 
driven gain {dM a /dt) of 35 Gta -1 for 2003-07, compared 
with a gain of 10 Gta -1 during 1992-2002 (Fig. 11). The 
derived ablation-dynamic loss rate (dA4 bd /df) of 206 Gta -1 
for 2003-07 is much larger than the corresponding loss rate 
of 17 Gta -1 during 1992-2002. Although both components 
of dA4/dfare larger in the later period, the increase in the loss 
term (-189 Gta -1 ) from melting and ice dynamics strongly 
dominates the increase in the gain term (+25 Gta -1 ) from 
snow accumulation. For comparison, if none of the mass 
change were ascribed to changes in accumulation, then net 
rate of mass loss for 2003-07 would be <1 71 Gta -1 (i.e. less 
negative than -171 Gta -1 ), because the positive dH a cA/dt 
component of observed elevation change would have an 
associated density of 0.91 instead of 0.527. 

To calculate a sensitivity of the derived dA 4/dt to our 
parameterization of A(t) f we ran the firn compaction for the 
2003-07 period using a parameterization of 8%5AK -1 
instead of 5%5AK -1 , which changed dA 4/dt by 
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Fig. 12. Components of mass change by drainage system. dM a /dt, dM bd /dt, d/Wdf (Gta -1 ) averaged over 500 m elevation bands for the eight 
drainage systems for 1992-2002 (black) and 2003-07 (red) with totals for 1992-2002 (black symbols) and 2003-07 (red symbols). 
Accumulation-driven mass increases are largest in DS3, DS4, DS7 and DS8, and dynamic-/ablation-driven thinning is largest in DS3, DS4 
and DS8 and very small in DS1, DS2, DS5 and DS6. 


11 Gta -1 , giving a sensitivity of -3.7 Gta -1 (%5A K) -1 . For 
2003-07, the sensitivity can also be estimated by noting that 
a dM a /dt value of 35 Gta -1 would be ~60Gta -1 (i.e. 
(35 x 0.91 )/0. 527) if part of the measured elevation changes 
were not ascribed to increasing A(t), giving a sensitivity of 
(35 - 60)/5 = -5.0 Gt a" 1 (% 5A K) -1 . For 1 992-2002, the dA4 a / 
dt of 10 Gta -1 would be 17 Gta -1 , giving -1.4 Gta -1 (%5 A 
K) -1 . Therefore, the sensitivity appears to depend somewhat 
on the specific A(T(t)) variability history. 

The relative values of the dA4 a /df and dM bd /dt com- 
ponents of dM/dt vary significantly from one drainage 
system to another, as do their changes with time (Fig. 12). 
In the northern drainage systems (DS1 and DS2), dM bd /dt 
values are essentially the same in both time periods, 
indicating little change in either the dynamics or the 
ablation in those systems. Also, dM bd /dt is not large 
compared with dA4 a /df in these drainage systems as it is in 
the southeast DS3 and DS4 and in the upper western DS7 
and DS8 (note the dM bd /dt scale is ten times the dA4 a /df scale 
in Fig. 12). Although dA4 a /df increases with time in all 
drainage systems, consistent with the increase in A(t), the 
increase in dM bd /dt from dynamic thinning and ablation 
dominates those increases in dM a /dt in all drainage systems 
except the northern DS1 and DS2. 

6. DISCUSSION 
6.1. Ice-sheet interior 

Our results show that the part of the accumulation zone 
above 2000 m (60.8% by area) is thickening on average 


during both periods by similar amounts, as shown by a dl/dt 
value of +4.78 cm a -1 during 2003-07 and the slightly larger 
value of +5.1 0cm a -1 during 1 992-2002 (Fig. 10). Although 
dl/dt\s similar in both periods, the accumulation-driven part 
(d/-/ a o\/df) increased from +1.13 to +4.68 cm a -1 while the 
dynamic driven part {dH d /dt) decreased from +3.96 to 
+0.1 0 cm a -1 . Therefore, the changes at higher elevations 
(>2000 m) are a combination of ice thickening of 
3.55 cm a -1 from increased A(t) and an approximately equal 
increase in dynamic thinning of 3.86 cm a -1 . 

If the resulting dl/dt decrease of 0.32 cm a -1 were not 
partitioned between dH a cA/dt and dhl d /dt, it would imply a 
decrease in the rate of mass gain of only 2 or 3 Gta -1 . 
However, the partitioning gives a dM a /dt increase of 
20 Gta -1 (i.e. gaining lower-density firn) and a larger dM d / 
dt decrease of 36 Gta -1 (i.e. losing higher-density ice), 
giving a net change of -16.5 Gta -1 above 2000 m. There- 
fore, although the accumulation-driven thickening approxi- 
mately balances the dynamic thinning (in terms of cm a -1 ), 
the net effect is a decrease in the rate of mass gain from 
44 to 28 Gta -1 . 

The map of the dM/dt distribution (Fig. 8), Table 1 and 
Figure 13 show that the northern DS1 and DS2 have net 
mass gains and have become slightly more positive during 
2003-07 compared with 1992-2002. The four drainage 
systems in the north and southwest (DS1, DS2, DS5 and 
DS6) show increased mass-gain rates above 2000 m (+18 to 
+31 Gta -1 ), suggesting increases in rates of ice accumu- 
lation, and small increased loss rates at lower elevations 
(-15 to -25 Gta -1 ). Overall, the major changes toward 
increased losses are in southeast DS3 and DS4 and west 


98 


Zwally and others: Greenland ice sheet mass balance 



Drainage subsystems 


Fig. 13. Total rate of mass change and rates above and below EL by 
drainage subsystems. 


central DS7 and DS8 in regions with accelerating glaciers 
(Rignot and Kanagaratnam, 2006). These four drainage 
systems have small negative changes above 2000m (+26 
to -2 Gta -1 ), but large changes at lower elevations (-32 to 
-1 56 Gta -1 ). 

6.2. Ice-sheet margins and outlet glaciers 

Only 1 1 .2% of the GIS area lies below the EL, and 39.2% 
lies below 2000 m, but the change in mass loss at these 
lower elevations, from -33 to -125 Gta -1 below the EL and 
from -51 to -198 Gta -1 below 2000 m (Table 1), dominates 
the total mass balance for the whole ice sheet. Therefore, the 
lower elevations of the GIS have been mostly thinning, and 
this thinning has increased significantly in the later period 
(Figs 10-12). EHowever, the increases below the EL in the 
northern DS1 and DS2 are small or near zero (Fig. 13). 

The only locations around the margins that show thicken- 
ing at lower elevations (Fig. 8) are the southeast corners of 
DS2.1 and DS2.2, a small area near Petermann Gletscher at 
the boundary between DS1.1 and DS1.2 and the small ice 
dome extending westward from the GIS at 66° N in DS6.2. In 
contrast, during 1992-2002, the pattern of thickening and 
thinning around the margin was more mixed, with areas of 
significant thickening in DS3.1 and DS5. Also, the thinning in 
the western ablation zone in DS6.2, DS7.2 and DS8.1 was 
small in 1992-2002 compared with 2003-07. 

About 25 Gta -1 or 5% of the total mass flux from the GIS is 
through Jakobshavn Isbrae (Rignot and Kanagaratnam, 2006), 
in DS7.1 . During the 1 990s, the mass gain of 5.8 Gta -1 above 
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Fig. 14. (a) d H/dt map of Jakobshavn Isbrae kriged on a 10 km grid 
from the along-track dh/dt plotted on the tracks, (b) d H/dt profile 
(interpolated from dH/dt map) and H profile along a central 
flowline (red curve in (a)). Red symbols are dh/dt± cr d/7 / df at points 
where the flowline crosses the ICESat tracks. Steep increase in dH/ 
df below the EL indicates effect of increased melting as well as 
thinning from glacier acceleration. 


the EL in DS7.1 exceeded the loss of 2.0 Gta -1 below the EL. 
For 2003-07, the gain above the EL reduced to 2.0 Gta -1 , 
and the loss below increased to 9.7 Gta -1 . Therefore, the net 
balance changed from a gain of 3.8 Gta -1 in the 1990s to a 
loss of 7.8 Gta -1 , which accounts for ~7% of the total 
increase in dAd/df between the two periods. 

Figures 14a, 15a, 16a and 17a show dH/dt maps of four 
outlet glaciers kriged from the along-track dh/dt at higher 
resolution (10 km, 5 km, 500 m and 500 m, respectively) 
than the 50 km used for dM/dt calculations and analysis. 
Figures 1 4b, 1 5b, 1 6b and 1 7b show a comparison between 
the dh/dt at intersections of a central flowline with the 
ICESat tracks and dH/dt interpolated from the higher- 
resolution maps. These maps provide a high-resolution 
pictorial representation of the topography on these outlet 
glaciers but are not used in our dM/dt calculations. 

The acceleration of Jakobshavn Isbrae (Joughin and others, 
2008a) as the floating ice tongue thinned and broke away 
during 1996-2004 certainly contributes to this increased 
rate of mass loss. However, the near-linear change in h(t) 
over 4 years, as shown in Figure 3c at six locations on 
Jakobshavn, suggests that the rate of thinning changed little 
during 2003-07. The rate of glacier thinning along the 
flowline (Fig. 14) shows a marked break at the EL (1440 m), 
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Fig. 15. (a) dH/dt map of Helheimgletscher kriged on a 5 km grid 
from along-track dh/dt plotted on the tracks, (b) dH/dt profile 
(interpolated from dH/dt map) and H profile along a central 
flowline (black curve in (a)). Black symbols are dh/dtdza dh/dt at 
points where the flowline crosses the ICESat tracks. Red and green 
segments in (a) indicate the ends of flowline. 


with increasingly larger rates of thinning at lower elevations 
downstream of the EL (slope of 1 .08 m a -1 (1 00 m) -1 . Above 
the EL the slope of the thinning rate is 0.38 m a -1 (1 00 m) -1 , 
or only one-third of that below the EL. Because the location 
of the EL is determined by the surface mass balance, and not 
by ice dynamics, this change at the EL cannot be explained 
by a change in the rate of dynamic thinning at this location. 
Assuming that the thinning above the EL is all caused by ice 
dynamics, extension of the 0.38 m a -1 (1 00 m) _1 thinning rate 
to below the EL suggests that at 1 000 m elevation, for 
example, about half of the 4ma _1 thinning is from 
dynamics and the other half is from increased melting. 
Therefore, increased melting and flow acceleration appear 
to be contributing similar amounts to the current thinning on 
the lower part of the glacier. Similarly, EHelheimgletscher on 
the east coast also shows a break in the thinning rate near 
the EL (Fig. 15). 

Two northern glaciers with floating ice tongues are 
Petermann Gletscher in the northeast corner of DS1.1 
(Fig. 16) and EHagen Brae in the northeast corner of DS1.3 
(Fig. 17). dH/dt on the floating part of Petermann shows an 
oscillating mixture of thickening and thinning, with an 
average close to zero, thinning inland of the grounding line 
to ~1 400 m, and thickening at higher elevations (Figs 1 6 and 
4c). The pattern of dH/dt on Hagen Brae is thickening on the 
floating part, strong thinning to about 600 m, and thickening 
inland to ~1 000 m. This pattern of thinning of grounded ice 
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Fig. 16. (a) dH/dt map of Petermann Gletscher kriged on a 500 m 
grid from along-track dh/dt plotted on the tracks, (b) dH/dt profile 
(interpolated from dH/dt map) and H profile along a central 
flowline (black curve in (a)). Black symbols are dh/dt± cr^h/dt at 
points where the flowline crosses the ICESat tracks. Red and green 
segments in (a) indicate the ends of flowline. 


inland of ice tongues that are not significantly thinning 
suggests that factors other than a decrease in the buttressing 
effect (Thomas and others, 2004) of floating 
ice are affecting the thinning of the grounded ice on 
these glaciers. 

6.3. Dynamic thinning and inland coupling 

The increase in d/-/ bd /df in 2003-07 compared with 1992- 
2002 (Fig. 10) is large below the EL (-19 to -72 cm a -1 ), 
where it is a combination of local ablation changes and 
dynamic thinning (Pritchard and others, 2009). Above the 
EL, the dynamic thinning (d/-/ d /df) decreases with increasing 
elevation and becomes zero or slightly positive above 
~2300m in averages over the GIS (Fig. lOe). The marked 
increase in dynamic thinning in 2003-07 extends into the 
GIS interior, decreasing with elevation (-6 to -30 cm a -1 
between 1 200 and 1 700 m and 5 to 1 cm a -1 above 2700 m). 

The increase in dynamic thinning during the later period 
indicates that the effects of increased melting and accel- 
erating glaciers, which are occurring at the margins, are 
dynamically extending inland to the highest elevations in 
some drainage systems. The increased thinning rate between 
1500 and 2000 m, for example, cannot be significantly 
affected dynamically by local changes in the surface mass 
balance, because the resulting changes in driving stresses 
from A(t) are small. Also, this inland extension occurs 
relatively fast (decadal-scale response). 
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Fig. 17 . (a) d H/dt map of Hagen Brae kriged on a 500 m grid from 
along-track dh/dt plotted on the tracks, (b) dH/dt profile (inter- 
polated from dH/dt map) and H profile along a central flowline 
(black curve in (a)). Black symbols are dh/dt ± cr^h/dt at points where 
the flowline crosses the ICESat tracks. Red and green segments in (a) 
indicate the ends of flowline. 

However, the inland extension of dynamic thinning is not 
the same in all drainage systems. The largest changes and 
inland extension are in DS4 in the southeast, followed by 
DS7, DS3, DS5, DS6 and DS8 (as shown by d/VIbd /df in 
Fig. 12). The response time for inland dynamic coupling 
should be fastest in DS4 (as observed; Fig. 12) where the ice 
is warmer and deforms more easily and where the slopes are 
steep and the distance from the margin to the ice divide is 
smaller, i.e. 100-200 km. In contrast, dH^/dt is nearly the 
same in both time periods in DS1 and DS2 in the north, 
where the ice is colder and the slopes are less steep. 

6.4. Response to climate change 

Changes in A(t) and T(t) on the GIS have immediate effects 
on the surface mass balance, the surface elevation and ice 
thickness, D, and the surface slope, a. Changes in D and a 
also change the driving stresses that force the ice flow, but 
these changes are generally expected to affect the ice 
dynamics only on timescales of centuries (Huybrechts and 
de Wolde, 1 999). Two mechanisms for rapid response of the 
ice dynamics to climate that have been observed in 
Greenland are changes in buttressing of outlet glaciers 
caused by removal of their floating ice tongues (Thomas and 
others, 2004; Nick and others, 2009) and changes in basal 
lubrication and sliding from the melt/acceleration effect 
(Zwally and others, 2002b; Joughin and others, 2008b; Van 
de Wal and others, 2008; Bartholomew and others, 2010). 


The observed increase of 1 64 Gt a -1 in d/Wdfto 1 71 Gta -1 
suggests that the GIS is responding to enhanced climate 
warming in the Arctic (AMAP, 2009), and that mass- 
exchange processes and ice dynamics causing greater mass 
loss are dominating those of mass gain. The thinning at the 
margins and thickening inland during the 1 990s was already 
consistent with a response to a warmer climate, even though 
the mass balance then was close to zero. From 1 997 to 2007, 
the 10 year trend in our AVHRR temperature record is 
+1 .7 K(10a) _1 , which compares with a 2.7 K(10a) -1 increase 
(Hall and others, 2008) for 2000-06 and a 1.8 K (10 a) -1 
increase (Box and others, 2009) for 1994-2007. Our corres- 
ponding increase in A(T(t)) is 8.5% (10 a) -1 , which compares 
with 7% (10 a)" 1 for 1988-2004 and 3.3% (10 a) -1 for 1958- 
2003 from regional-climate surface mass-balance models 
(Hanna and others, 2005; Box and others, 2006). 

The increased mass loss is partly due to changes in 
surface mass balance and partly to changes in the ice 
dynamics. Since temperature increase is the primary forcing 
for these changes, including the variation of A(t) with T(t), it 
might be expected that an additional warming of 3.4 K in 
another two decades would increase the rate of loss to 
500 Gta -1 (1.4 mm a -1 sea-level rise), assuming a linear 
relationship between temperature increase and mass loss. 
Numerical models of ice-sheet response to climate change 
in the next century generally account for changes in surface 
mass balance, but either neglect the effects of surface 
changes on the ice dynamics or only consider local changes 
in the driving forces. Most models do not include higher- 
order stress couplings within the ice, melt/acceleration 
effects, or the acceleration of outlet glaciers. The observed 
increase in dynamic thinning in 2003-07 and the dynamic 
coupling inland on decadal timescales is a new indication of 
how the GIS is responding to climate change that needs to 
be accounted for in predicting ice-sheet behavior. 

7. SUMMARY 

The mass changes of the GIS that we derive for two periods 
(1992-2002 and 2003-07) show a significant increase in 
the rate of mass loss (from 7±3Gta -1 to 171 ±4 Gta -1 ) 
during a period of significant climate warming in Green- 
land (~2K(10a) -1 ). From 1992-2002 to 2003-07, the 
contribution of the GIS to sea-level rise changed from 
essentially zero to 0.5 mm a -1 , ~15% of the recent rate of 
sea-level rise (Solomon and others, 2007). Even though the 
GIS was close to mass balance during the 1990s, it was 
already showing characteristics of responding to a warmer 
climate, specifically thinning at the margins and thickening 
inland at higher elevations. Our results show that increased 
ice thinning due to increases in melting and acceleration of 
outlet glaciers during 2003-07 now strongly exceeds the 
inland thickening from increases in accumulation. Over the 
entire GIS, the mass loss between the two periods, from 
increased melting and ice dynamics, increased from 17 to 
206 Gta -1 while the mass gain, from increased precipi- 
tation and accumulation, increased from 10 to 25 Gta -1 . 
These results provide new quantitative information on how 
the competing mass-exchange processes that affect ice- 
sheet growth or shrinkage evolve with time in a changing 
climate. Understanding this evolution is important for the 
future, because the climate affecting the Greenland and 
Antarctic ice sheets is predicted to continue warming as the 
build-up of C0 2 and other greenhouse gases in the Earth's 
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atmosphere continues at an increasing rate (Solomon and 
others, 2007). 

Our separation of the mass changes, dM a /dt, driven by 
accumulation variability, from the changes, dA4>d /dt, driven 
by changes in ablation and ice dynamics is enabled by a 
new formulation of the elevation changes as perturbations 
from steady state. This formulation, as applied to measure- 
ments over two time periods (1 992-2002 and 2003-07) with 
midpoints separated by 8 years, also provides a new 
determination of the changes in the dynamic thickening or 
dynamic thinning, and how the changes at the margins are 
propagating inland on decadal timescales. 

Our detailed distributions derived for dH a cA/dt and dH^/ 
dt in Figure 7a and b and for dM a /dt, dM bd /dt and dM/dt (as a 
function of elevation and by drainage system in Figure 12) 
show how these parameters vary significantly among the 
different climate regimes of the GIS. The eight drainage 
systems have significant differences in important parameters 
that affect the surface balance and ice dynamics, including 
surface temperature, accumulation rate, internal ice tem- 
perature, basal melting or freezing, number of outlet glaciers 
with or without floating ice tongues, and the length of the ice 
margin terminating on land versus ocean. For example, the 
northern DS1 and DS2 have lower temperatures, lower 
accumulation rates and smaller ablation zones than the more 
southerly systems. Also, the drainage systems on the west 
coast have wide ablation zones with low slopes, and most of 
DS5 and DS6 terminate on land, whereas the east central 
DS3, the southeast DS4 and the west central DS7 and DS8 
have many outlet glaciers terminating in the ocean and some 
with floating ice tongues. Our results provide a new 
description of the behavior of the GIS in these differing 
climate and glaciological regimes and therefore an improved 
basis for testing models of the ice response to climate change. 
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